home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 2 / Apprentice-Release2.iso / Source Code / C / Utilities / Calc 1.24.7 / lib / lucas.cal < prev    next >
Encoding:
Text File  |  1992-02-24  |  27.2 KB  |  1,015 lines  |  [TEXT/R*ch]
  1. /*
  2.  * Copyright (c) 1992 Landon Curt Noll
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * By: Landon Curt Noll
  7.  *     chongo@toad.com  -or-  ...!{pyramid,sun,uunet}!sun!hoptoad!chongo
  8.  *
  9.  *
  10.  * lucas - perform a Lucas primality test on h*2^n-1
  11.  *
  12.  * HISTORICAL NOTE:
  13.  *
  14.  * On 6 August 1989 at 00:53 PDT, the 'Amdahl 6', a team consisting of
  15.  * John Brown, Landon Curt Noll, Bodo Parady, Gene Smith, Joel Smith and
  16.  * Sergio Zarantonello proved the following 65087 digit number to be prime:
  17.  *
  18.  *                  216193
  19.  *            391581 * 2      -1
  20.  *
  21.  * The primality was demonstrated by a program implementing the test
  22.  * found in these routines.  An Amdahl 1200 takes 1987 seconds to test
  23.  * the primality of this number.  A Cray 2 took several hours to
  24.  * confirm this prime.  As of 23 Jul 1991, this prime was still the
  25.  * largest known prime.
  26.  *
  27.  * The same team also discovered the following twin prime pair:
  28.  *
  29.  *               11235           11235
  30.  *        1706595 * 2     -1    1706595 * 2     +1
  31.  *
  32.  * As of 23 Jul 1991, these primes was still the largest known twin prime pair.
  33.  *
  34.  * ON GAINING A WORLD RECORD:
  35.  *
  36.  * The routines in calc were designed to be portable, and to work on
  37.  * numbers of 'sane' size.  The 'ultra-high speed multi-precision'
  38.  * package was a highly machine dependent collection of routines tuned
  39.  * to work with very large numbers.  The heart of this package was a
  40.  * multiplication and square routine that were based on Fast Fourier
  41.  * Transforms.  Details of the FFT are to be published in a up-coming
  42.  * paper to be written by the 'Amdahl 6'.
  43.  *
  44.  * Having a fast computer, and a good multi-precision package are
  45.  * critical, but one also needs to know where to look in order to have
  46.  * a good chance at a record.  Knowing what to test is beyond the scope
  47.  * of this routine.  However the following observations are noted:
  48.  *
  49.  *    test numbers of the form h*2^n-1
  50.  *    fix a value of n and vary the value h
  51.  *    n mod 128 == 0
  52.  *    h*2^n-1 is not divisible by any small prime < 2^40
  53.  *    0 < h < 2^39
  54.  *    h*2^n+1 is not divisible by any small prime < 2^40
  55.  *
  56.  * The Mersenne test for '2^n-1' is the fastest known primality test
  57.  * for a given large numbers.  However, it is faster to search for
  58.  * primes of the form 'h*2^n-1'.  When n is around 20000, one can find
  59.  * a prime of the form 'h*2^n-1' in about 1/2 the time.
  60.  *
  61.  * Critical to understanding why 'h*2^n-1' is to observe that primes of
  62.  * the form '2^n-1' seem to bunch around "islands".  Such "islands"
  63.  * seem to be getting fewer and farther in-between, forcing the time
  64.  * for each test to grow longer and longer (worse then O(n^2 log n)).
  65.  * On the other hand, when one tests 'h*2^n-1', fixes 'n' and varies
  66.  * 'h', the time to test each number remains relatively constant.
  67.  *
  68.  * It is clearly a win to eliminate potential test candidates by
  69.  * rejecting numbers that that are divisible by 'small' primes.  We
  70.  * (the "Amdahl 6") rejected all numbers that were divisible by primes
  71.  * less than '2^40'.  We stopped looking for small factors at '2^40'
  72.  * when the rate of candidates being eliminated was slowed down to
  73.  * just a trickle.
  74.  *
  75.  * The 'n mod 128 == 0' restriction allows one to test for divisability
  76.  * of small primes more quickly.  To test of 'q' is a factor of 'k*2^n-1',
  77.  * one check to see if 'k*2^n mod q' == 1, which is the same a checking
  78.  * if 'h*(2^n mod q) mod q' == 1.  One can compute '2^n mod q' by making
  79.  * use of the following:
  80.  *
  81.  *    if
  82.  *        y = 2^x mod q
  83.  *    then
  84.  *        2^(2x) mod q   == y^2 mod q        0 bit
  85.  *        2^(2x+1) mod q == 2*y^2 mod q        1 bit
  86.  *
  87.  * The choice of which expression depends on the binary pattern of 'n'.
  88.  * Since '1' bits require an extra step (multiply by 2), one should
  89.  * select value of 'n' that contain mostly '0' bits.  The restriction
  90.  * of 'n mod 128 == 0' ensures that the bottom 7 bits of 'n' are 0.
  91.  *
  92.  * By limiting 'h' to '2^39' and eliminating all values divisible by
  93.  * small primes < twice the 'h' limit (2^40), one knows that all
  94.  * remaining candidates are relatively prime.  Thus, when a candidate
  95.  * is proven to be composite (not prime) by the big test, one knows
  96.  * that the factors for that number (whatever they may be) will not
  97.  * be the factors of another candidate.
  98.  *
  99.  * Finally, one should eliminate all values of 'h*2^n-1' where
  100.  * 'h*2^n+1' is divisable by a small primes.  The ideas behind this 
  101.  * point is beyond the scope of this program and will be discussed
  102.  * in the same up-comming paper.
  103.  */
  104.  
  105. global pprod256;    /* product of  "primes up to 256" / "primes up to 46" */
  106. pprod256 = 0;
  107.  
  108. dbg = 0;        /* 1 => print debug statements */
  109.  
  110. /*
  111.  * lucas - lucas primality test on h*2^n-1
  112.  *
  113.  * ABOUT THE TEST:
  114.  *
  115.  * This routine will perform a primality test on h*2^n-1 based on
  116.  * the mathematics of Lucas, Lehmer and Riesel.  One should read
  117.  * the following article:
  118.  *
  119.  * Ref1:
  120.  *    "Lucasian Criteria for the Primality of N=h*2^n-1", by Hans Riesel,
  121.  *    Mathematics of Computation, Vol 23 #108, pp. 869-875, Oct 1969
  122.  *
  123.  * The following book is also useful:
  124.  *
  125.  * Ref2:
  126.  *    "Prime numbers and Computer Methods for Factorization", by Hans Riesel,
  127.  *    Birkhauser, 1985, pp 131-134, 278-285, 438-444
  128.  *
  129.  * A few useful Legendre identities may be found in:
  130.  *
  131.  * Ref3:
  132.  *    "Introduction to Analytic Number Theory", by Tom A. Apostol,
  133.  *    Springer-Verlag, 1984, p 188.
  134.  *
  135.  * This test is performed as follows:    (see Ref1, Theorem 5)
  136.  *
  137.  *    a) generate u(0)         (see the function gen_u0() below)
  138.  *
  139.  *    b) generate u(n-2) according to the rule:
  140.  *
  141.  *        u(i+1) = u(i)^2-2 mod h*2^n-1
  142.  *
  143.  *    c) h*2^n-1 is prime if and only if u(n-2) == 0        Q.E.D. :-)
  144.  *
  145.  * Now the following conditions must be true for the test to work:
  146.  *
  147.  *      n >= 2
  148.  *    h >= 1
  149.  *      h < 2^n
  150.  *    h mod 2 == 1
  151.  *
  152.  * A few misc notes:
  153.  *
  154.  * In order to reduce the number of tests, as attempt to eliminate
  155.  * any number that is divisible by a prime less than 257.  Valid prime
  156.  * candidates less than 257 are declared prime as a special case.
  157.  *
  158.  * The condition 'h mod 2 == 1' is not a problem.  Say one is testing
  159.  * 'j*2^m-1', where j is even.  If we note that:
  160.  *
  161.  *      j mod 2^x == 0 for x>0   implies   j*2^m-1 == ((j/2^x)*2^(m+x))-1,
  162.  *
  163.  * then we can let h=j/2^x and n=m+x and test 'h*2^n-1' which is the value.
  164.  * We need only consider odd values of h because we can rewrite our numbers
  165.  * do make this so.
  166.  *
  167.  * input:
  168.  *    h    the h as in h*2^n-1
  169.  *    n    the n as in h*2^n-1
  170.  *
  171.  * returns:
  172.  *    1 => h*2^n-1 is prime
  173.  *    0 => h*2^n-1 is not prime
  174.  *     -1 => a test could not be formed, or h >= 2^n, h <= 0, n <= 0
  175.  */
  176. define
  177. lucas(h, n)
  178. {
  179.     local testval;        /* h*2^n-1 */
  180.     local shiftdown;    /* the power of 2 that divides h */
  181.     local u;        /* the u(i) sequence value */
  182.     local v1;        /* the v(1) generator of u(0) */
  183.     local i;        /* u sequence cycle number */
  184.     local oldh;        /* pre-reduced h */
  185.     local oldn;        /* pre-reduced n */
  186.     local bits;        /* highbit of h*2^n-1 */
  187.  
  188.     /*
  189.      * check arg types
  190.      */
  191.     if (!isint(h)) {
  192.         ldebug("lucas", "h is non-int");
  193.         quit "FATAL: bad args: h must be an integer";
  194.     }
  195.     if (!isint(n)) {
  196.         ldebug("lucas", "n is non-int");
  197.         quit "FATAL: bad args: n must be an integer";
  198.     }
  199.  
  200.     /*
  201.      * reduce h if even
  202.      *
  203.      * we will force h to be odd by moving powers of two over to 2^n
  204.      */
  205.     oldh = h;
  206.     oldn = n;
  207.     shiftdown = fcnt(h,2);  /* h % 2^shiftdown == 0, max shiftdown */
  208.     if (shiftdown > 0) {
  209.         h >>= shiftdown;
  210.         n += shiftdown;
  211.     }
  212.  
  213.     /*
  214.      * enforce the 0 < h < 2^n rule
  215.      */
  216.     if (h <= 0 || n <= 0) {
  217.         print "ERROR: reduced args violate the rule: 0 < h < 2^n";
  218.         print "    ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n;
  219.         ldebug("lucas", "unknown: h <= 0 || n <= 0");
  220.         return -1;
  221.     }
  222.     if (highbit(h) >= n) {
  223.         print "ERROR: reduced args violate the rule: h < 2^n";
  224.         print "    ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n;
  225.         ldebug("lucas", "unknown: highbit(h) >= n");
  226.         return -1;
  227.     }
  228.  
  229.     /*
  230.      * catch the degenerate case of h*2^n-1 == 1
  231.      */
  232.     if (h == 1 && n == 1) {
  233.         ldebug("lucas", "not prime: h == 1 && n == 1");
  234.         return 0;    /* 1*2^1-1 == 1 is not prime */
  235.     }
  236.  
  237.     /*
  238.      * catch the degenerate case of n==2
  239.      *
  240.      * n==2 and 0<h<2^n  ==>  0<h<4
  241.      *
  242.      * Since h is now odd  ==>  h==1 or h==3
  243.      */
  244.     if (h == 1 && n == 2) {
  245.         ldebug("lucas", "prime: h == 1 && n == 2");
  246.         return 1;    /* 1*2^2-1 == 3 is prime */
  247.     }
  248.     if (h == 3 && n == 2) {
  249.         ldebug("lucas", "prime: h == 3 && n == 2");
  250.         return 1;    /* 3*2^2-1 == 11 is prime */
  251.     }
  252.  
  253.     /*
  254.      * catch small primes < 257
  255.      *
  256.      * We check for only a few primes because the other primes < 257
  257.      * violate the checks above.
  258.      */
  259.     if (h == 1) {
  260.         if (n == 3 || n == 5 || n == 7) {
  261.             ldebug("lucas", "prime: 3, 7, 31, 127 are prime");
  262.             return 1;    /* 3, 7, 31, 127 are prime */
  263.         }
  264.     }
  265.     if (h == 3) {
  266.         if (n == 2 || n == 3 || n == 4 || n == 6) {
  267.             ldebug("lucas", "prime: 11, 23, 47, 191 are prime");
  268.             return 1;    /* 11, 23, 47, 191 are prime */
  269.         }
  270.     }
  271.     if (h == 5 && n == 4) {
  272.         ldebug("lucas", "prime: 79 is prime");
  273.         return 1;        /* 79 is prime */
  274.     }
  275.     if (h == 7 && n == 5) {
  276.         ldebug("lucas", "prime: 223 is prime");
  277.         return 1;        /* 223 is prime */
  278.     }
  279.     if (h == 15 && n == 4) {
  280.         ldebug("lucas", "prime: 239 is prime");
  281.         return 1;        /* 239 is prime */
  282.     }
  283.  
  284.     /*
  285.      * Avoid any numbers divisable by small primes
  286.      */
  287.     /*
  288.      * check for 3 <= prime factors < 29
  289.      * pfact(28)/2 = 111546435
  290.      */
  291.     testval = h*2^n - 1;
  292.     if (gcd(testval, 111546435) > 1) {
  293.         /* a small 3 <= prime < 29 divides h*2^n-1 */
  294.         ldebug("lucas","not-prime: 3<=prime<29 divides h*2^n-1");
  295.         return 0;
  296.     }
  297.     /*
  298.      * check for 29 <= prime factors < 47
  299.      * pfact(46)/pfact(28) = 5864229
  300.      */
  301.     if (gcd(testval, 58642669) > 1) {
  302.         /* a small 29 <= prime < 47 divides h*2^n-1 */
  303.         ldebug("lucas","not-prime: 29<=prime<47 divides h*2^n-1");
  304.         return 0;
  305.     }
  306.     /*
  307.      * check for prime 47 <= factors < 257, if h*2^n-1 is large
  308.      * 2^282 > pfact(256)/pfact(46) > 2^281
  309.      */
  310.     bits = highbit(testval);
  311.     if (bits >= 281) {
  312.         if (pprod256 <= 0) {
  313.             pprod256 = pfact(256)/pfact(46);
  314.         }
  315.         if (gcd(testval, pprod256) > 1) {
  316.             /* a small 47 <= prime < 257 divides h*2^n-1 */
  317.             ldebug("lucas",\
  318.                 "not-prime: 47<=prime<257 divides h*2^n-1");
  319.             return 0;
  320.         }
  321.     }
  322.  
  323.     /*
  324.      * try to compute u(0)
  325.      *
  326.      * We will use gen_v1() to give us a v(1) using the values
  327.      * of 'h' and 'n'.  We will then use gen_u0() to convert
  328.      * the v(1) into u(0).
  329.      *
  330.      * If gen_v1() returns a negative value, then we failed to
  331.      * generate a test for h*2^n-1.  This is because h mod 3 == 0
  332.      * is hard to do, and in rare cases, exceed the tables found
  333.      * in this program.  We will generate an message and assume
  334.      * the number is not prime, even though if we had a larger
  335.      * table, we might have been able to show that it is prime.
  336.      */
  337.     v1 = gen_v1(h, n, testval);
  338.     if (v1 < 0) {
  339.         /* failure to test number */
  340.         print "unable to compute v(1) for", h : "*2^" : n : "-1";
  341.         ldebug("lucas", "unknown: no v(1)");
  342.         return -1;
  343.     }
  344.     u = gen_u0(h, n, testval, v1);
  345.  
  346.     /*
  347.      * compute u(n-2)
  348.      */
  349.     for (i=3; i <= n; ++i) {
  350.         u = (u^2 - 2) % testval;
  351.     }
  352.  
  353.     /*
  354.      * return 1 if prime, 0 is not prime
  355.      */
  356.     if (u == 0) {
  357.         ldebug("lucas", "prime: end of test");
  358.         return 1;
  359.     } else {
  360.         ldebug("lucas", "not-prime: end of test");
  361.         return 0;
  362.     }
  363. }
  364.  
  365. /*
  366.  * gen_u0 - determine the initial Lucas sequence for h*2^n-1
  367.  *
  368.  * According to Ref1, Theorem 5:
  369.  *
  370.  *    u(0) = alpha^h + alpha^(-h)
  371.  *
  372.  * Now:
  373.  *
  374.  *    v(x) = alpha^x + alpha^(-x)    (Ref1, bottom of page 872)
  375.  *
  376.  * Therefore:
  377.  *
  378.  *    u(0) = v(h)
  379.  *
  380.  * We calculate v(h) as follows:    (Ref1, top of page 873)
  381.  *
  382.  *    v(0) = alpha^0 + alpha^(-0) = 2
  383.  *    v(1) = alpha^1 + alpha^(-1) = gen_v1(h,n)
  384.  *    v(n+2) = v(1)*v(n+1) - v(n)
  385.  *
  386.  * This function does not concern itself with the value of 'alpha'.
  387.  * The gen_v1() function is used to compute v(1), and identity
  388.  * functions take it from there.
  389.  *
  390.  * It can be shown that the following are true:
  391.  *
  392.  *    v(2*n) = v(n)^2 - 2
  393.  *    v(2*n+1) = v(n+1)*v(n) - v(1)
  394.  *
  395.  * To prevent v(x) from growing too large, one may replace v(x) with
  396.  * `v(x) mod h*2^n-1' at any time.
  397.  *
  398.  * See the function gen_v1() for details on the value of v(1).
  399.  *
  400.  * input:
  401.  *    h    - h as in h*2^n-1    (h mod 2 != 0)
  402.  *    n    - n as in h*2^n-1
  403.  *    testval    - h*2^n-1
  404.  *    v1    - gen_v1(h,n)        (see function below)
  405.  *
  406.  * returns:
  407.  *    u(0)    - initial value for Lucas test on h*2^n-1
  408.  *    -1    - failed to generate u(0)
  409.  */
  410. define
  411. gen_u0(h, n, testval, v1)
  412. {
  413.     local shiftdown;    /* the power of 2 that divides h */
  414.     local r;        /* low value: v(n) */
  415.     local s;        /* high value: v(n+1) */
  416.     local hbits;        /* highest bit set in h */
  417.     local i;
  418.  
  419.     /*
  420.      * check arg types
  421.      */
  422.     if (!isint(h)) {
  423.         quit "bad args: h must be an integer";
  424.     }
  425.     if (!isint(n)) {
  426.         quit "bad args: n must be an integer";
  427.     }
  428.     if (!isint(testval)) {
  429.         quit "bad args: testval must be an integer";
  430.     }
  431.     if (!isint(v1)) {
  432.         quit "bad args: v1 must be an integer";
  433.     }
  434.     if (testval <= 0) {
  435.         quit "bogus arg: testval is <= 0";
  436.     }
  437.     if (v1 <= 0) {
  438.         quit "bogus arg: v1 is <= 0";
  439.     }
  440.  
  441.     /*
  442.      * enforce the h mod rules
  443.      */
  444.     if (h%2 == 0) {
  445.         quit "h must not be even";
  446.     }
  447.  
  448.     /*
  449.      * enforce the h > 0 and n >= 2 rules
  450.      */
  451.     if (h <= 0 || n < 1) {
  452.         quit "reduced args violate the rule: 0 < h < 2^n";
  453.     }
  454.     hbits = highbit(h);
  455.     if (hbits >= n) {
  456.         quit "reduced args violate the rule: 0 < h < 2^n";
  457.     }
  458.  
  459.     /*
  460.      * build up u2 based on the reversed bits of h
  461.      */
  462.     /* setup for bit loop */
  463.     r = v1;
  464.     s = (r^2 - 2);
  465.  
  466.     /*
  467.      * deal with small h as a special case
  468.      *
  469.      * The h value is odd > 0, and it needs to be
  470.      * at least 2 bits long for the loop below to work.
  471.      */
  472.     if (h == 1) {
  473.         ldebug("gen_u0", "quick h == 1 case");
  474.         return r%testval;
  475.     }
  476.  
  477.     /* cycle from second highest bit to second lowest bit of h */
  478.     for (i=hbits-1; i > 0; --i) {
  479.  
  480.         /* bit(i) is 1 */
  481.         if (isset(h,i)) {
  482.  
  483.             /* compute v(2n+1) = v(r+1)*v(r)-v1 */
  484.             r = (r*s - v1) % testval;
  485.  
  486.             /* compute v(2n+2) = v(r+1)^2-2 */
  487.             s = (s^2 - 2) % testval;
  488.  
  489.         /* bit(i) is 0 */
  490.         } else {
  491.  
  492.             /* compute v(2n+1) = v(r+1)*v(r)-v1 */
  493.             s = (r*s - v1) % testval;
  494.  
  495.             /* compute v(2n) = v(r)^-2 */
  496.             r = (r^2 - 2) % testval;
  497.         }
  498.     }
  499.  
  500.     /* we know that h is odd, so the final bit(0) is 1 */
  501.     r = (r*s - v1) % testval;
  502.  
  503.     /* compute the final u2 return value */
  504.     return r;
  505. }
  506.  
  507. /*
  508.  * Trial tables used by gen_v1()
  509.  *
  510.  * When h mod 3 == 0, one needs particular values of D, a and b (see gen_v1
  511.  * documentation) in order to find a value of v(1).
  512.  *
  513.  * This table defines 'quickmax' possible tests to be taken in ascending
  514.  * order.  The v1_qval[x] refers to a v(1) value from Ref1, Table 1.  A
  515.  * related D value is found in d_qval[x].  All D values expect d_qval[1]
  516.  * are also taken from Ref1, Table 1.  The case of D == 21 as listed in
  517.  * Ref1, Table 1 can be changed to D == 7 for the sake of the test because
  518.  * of {note 6}.
  519.  *
  520.  * It should be noted that the D values all satisfy the selection values
  521.  * as outlined in the gen_v1() function comments.  That is:
  522.  *
  523.  *       D == P*(2^f)*(3^g)
  524.  *
  525.  * where f == 0 and g == 0, P == D.  So we simply need to check that
  526.  * one of the following two cases are true:
  527.  *
  528.  *       P mod 4 ==  1  and  J(h*2^n-1 mod P, P) == -1
  529.  *       P mod 4 == -1  and  J(h*2^n-1 mod P, P) ==  1
  530.  *
  531.  * In all cases expect one case, the value of r is:
  532.  *
  533.  *       r == Q*(2^j)*(3^k)*(z^2)
  534.  *
  535.  * where Q == 1.  In these cases, no further checking is needed,
  536.  * and the v(1) value can be returned.
  537.  */
  538. quickmax = 8;
  539. mat d_qval[quickmax];
  540. mat v1_qval[quickmax];
  541. d_qval[0] = 5;        v1_qval[0] = 3;        /* a=1   b=1  r=4  */
  542. d_qval[1] = 7;        v1_qval[1] = 5;        /* a=3   b=1  r=12  D=21 */
  543. d_qval[2] = 13;        v1_qval[2] = 11;    /* a=3   b=1  r=4  */
  544. d_qval[3] = 11;        v1_qval[3] = 20;    /* a=3   b=1  r=2  */
  545. d_qval[4] = 29;        v1_qval[4] = 27;    /* a=5   b=1  r=4  */
  546. d_qval[5] = 53;        v1_qval[5] = 51;    /* a=53  b=1  r=4  */
  547. d_qval[6] = 17;        v1_qval[6] = 66;    /* a=17  b=1  r=1  */
  548. d_qval[7] = 19;        v1_qval[7] = 74;    /* a=38  b=1  r=2  */
  549.  
  550. /*
  551.  * gen_v1 - compute the v(1) for a given h*2^n-1 if we can
  552.  *
  553.  * This function assumes:
  554.  *
  555.  *    n > 2            (n==2 has already been eliminated)
  556.  *    h mod 2 == 1
  557.  *    h < 2^n
  558.  *    h*2^n-1 mod 3 != 0    (h*2^n-1 has no small factors, such as 3)
  559.  *
  560.  * The generation of v(1) depends on the value of h.  There are two cases
  561.  * to consider, h mod 3 != 0, and h mod 3 == 0.
  562.  *
  563.  ***
  564.  *
  565.  * Case 1:    (h mod 3 != 0)
  566.  *
  567.  * This case is easy and always finds v(1).
  568.  *
  569.  * In Ref1, page 869, one finds that if:    (or see Ref2, page 131-132)
  570.  *
  571.  *    h mod 6 == +/-1
  572.  *    h*2^n-1 mod 3 != 0
  573.  *
  574.  * which translates, gives the functions assumptions, into the condition:
  575.  *
  576.  *    h mod 3 != 0
  577.  *
  578.  * If this case condition is true, then:
  579.  *
  580.  *    u(0) = (2+sqrt(3))^h + (2-sqrt(3))^h        (see Ref1, page 869)
  581.  *         = (2+sqrt(3))^h + (2+sqrt(3))^(-h)
  582.  *
  583.  * and since Ref1, Theorem 5 states:
  584.  *
  585.  *    u(0) = alpha^h + alpha^(-h)
  586.  *    r = abs(2^2 - 1^2*3) = 1
  587.  *
  588.  * and the bottom of Ref1, page 872 states:
  589.  *
  590.  *    v(x) = alpha^x + alpha^(-x)
  591.  *
  592.  * If we let:
  593.  *
  594.  *    alpha = (2+sqrt(3))
  595.  *
  596.  * then
  597.  *
  598.  *    u(0) = v(h)
  599.  *
  600.  * so we simply return
  601.  *
  602.  *    v(1) = alpha^1 + alpha^(-1)
  603.  *         = (2+sqrt(3)) + (2-sqrt(3))
  604.  *         = 4
  605.  *
  606.  ***
  607.  *
  608.  * Case 2:    (h mod 3 == 0)
  609.  *
  610.  * This case is not so easy and finds v(1) in most all cases.  In this
  611.  * version of this program, we will simply return -1 (failure) if we
  612.  * hit one of the cases that fall thru the cracks.  This does not happen
  613.  * often, so this is not too bad.
  614.  *
  615.  * Ref1, Theorem 5 contains the following definitions:
  616.  *
  617.  *        r = abs(a^2 - b^2*D)
  618.  *        alpha = (a + b*sqrt(D))^2/r
  619.  *
  620.  * where D is 'square free', and 'alpha = epsilon^s' (for some s>0) are units
  621.  * in the quadratic field K(sqrt(D)).
  622.  *
  623.  * One can find possible values for a, b and D in Ref1, Table 1 (page 872).
  624.  * (see the file lucas_tbl.cal)
  625.  *
  626.  * Now Ref1, Theorem 5 states that if:
  627.  *
  628.  *    L(D, h*2^n-1) = -1                [condition 1]
  629.  *    L(r, h*2^n-1) * (a^2 - b^2*D)/r = -1        [condition 2]
  630.  *
  631.  * where L(x,y) is the Legendre symbol (see below), then:
  632.  *
  633.  *    u(0) = alpha^h + alpha^(-h)
  634.  *
  635.  * The bottom of Ref1, page 872 states:
  636.  *
  637.  *    v(x) = alpha^x + alpha^(-x)
  638.  *
  639.  * thus since:
  640.  *
  641.  *    u(0) = v(h)
  642.  *
  643.  * so we want to return:
  644.  *
  645.  *    v(1) = alpha^1 + alpha^(-1)
  646.  *
  647.  * Therefore we need to take a given (D,a,b), determine if the two conditions
  648.  * are true, and return the related v(1).
  649.  *
  650.  * Before we address the two conditions, we need some background information
  651.  * on two symbols, Legendre and Jacobi.  In Ref 2, pp 278, 284-285, we find
  652.  * the following definitions of J(a,p) and L(a,n):
  653.  *
  654.  * The Legendre symbol L(a,p) takes the value:
  655.  *
  656.  *    L(a,p) == 1    => a is a quadratic residue of p
  657.  *    L(a,p) == -1    => a is NOT a quadratic residue of p
  658.  *
  659.  * when
  660.  *
  661.  *    p is prime
  662.  *    p mod 2 == 1
  663.  *    gcd(a,p) == 1
  664.  *
  665.  * The value x is a quadratic residue of y if there exists some integer z
  666.  * such that:
  667.  *
  668.  *    z^2 mod y == x
  669.  *
  670.  * The Jacobi symbol J(x,y) takes the value:
  671.  *
  672.  *    J(x,y) == 1    => y is not prime, or x is a quadratic residue of y
  673.  *    J(x,y) == -1    => x is NOT a quadratic residue of y
  674.  *
  675.  * when
  676.  *
  677.  *    y mod 2 == 1
  678.  *    gcd(x,y) == 1
  679.  *
  680.  * In the following comments on Legendre and Jacobi identities, we shall
  681.  * assume that the arguments to the symbolic are valid over the symbol
  682.  * definitions as stated above.
  683.  *
  684.  * In Ref2, pp 280-284, we find that:
  685.  *
  686.  *    L(a,p)*L(b,p) == L(a*b,p)                {A3.5}
  687.  *    J(x,y)*J(z,y) == J(x*z,y)                {A3.14}
  688.  *    L(a,p) == L(p,a) * (-1)^((a-1)*(p-1)/4)            {A3.8}
  689.  *    J(x,y) == J(y,x) * (-1)^((x-1)*(y-1)/4)            {A3.17}
  690.  *
  691.  * The equality L(a,p) == J(a,p) when:                {note 0}
  692.  *
  693.  *    p is prime
  694.  *    p mod 2 == 1
  695.  *    gcd(a,p) == 1
  696.  *
  697.  * It can be shown that (see Ref3):
  698.  *
  699.  *    L(a,p) == L(a mod p, p)                    {note 1}
  700.  *    L(z^2, p) == 1                        {note 2}
  701.  *
  702.  * From Ref2, table 32:
  703.  *
  704.  *    p mod 8 == +/-1   implies  L(2,p) == 1            {note 3}
  705.  *    p mod 12 == +/-1  implies  L(3,p) == 1            {note 4}
  706.  *
  707.  * Since h*2^n-1 mod 8 == -1, for n>2, note 3 implies:
  708.  *
  709.  *    L(2, h*2^n-1) == 1            (n>2)        {note 5}
  710.  *
  711.  * Since h=3*A, h*2^n-1 mod 12 == -1, for A>0, note 4 implies:
  712.  *
  713.  *    L(3, h*2^n-1) == 1                    {note 6}
  714.  *
  715.  * By use of {A3.5}, {note 2}, {note 5} and {note 6}, one can show:
  716.  *
  717.  *    L((2^g)*(3^l)*(z^2), h*2^n-1) == 1  (g>=0,l>=0,z>0,n>2)    {note 7}
  718.  *
  719.  * Returning to the testing of conditions, take condition 1:
  720.  *
  721.  *    L(D, h*2^n-1) == -1            [condition 1]
  722.  *
  723.  * In order for J(D, h*2^n-1) to be defined, we must ensure that D
  724.  * is not a factor of h*2^n-1.  This is done by pre-screening h*2^n-1 to
  725.  * not have small factors and selecting D less than that factor check limit.
  726.  *
  727.  * By use of {note 7}, we can show that when we choose D to be:
  728.  *
  729.  *    D is square free
  730.  *    D = P*(2^f)*(3^g)            (P is prime>2)
  731.  *
  732.  * The square free condition implies f = 0 or 1, g = 0 or 1.  If f and g
  733.  * are both 1, P must be a prime > 3.
  734.  *
  735.  * So given such a D value:
  736.  *
  737.  *    L(D, h*2^n-1) == L(P*(2^g)*(3^l), h*2^n-1)
  738.  *              == L(P, h*2^n-1) * L((2^g)*(3^l), h*2^n-1)       {A3.5}
  739.  *              == L(P, h*2^n-1) * 1                   {note 7}
  740.  *              == L(h*2^n-1, P)*(-1)^((h*2^n-2)*(P-1)/4)           {A3.8}
  741.  *              == L(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4)  {note 1}
  742.  *              == J(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4)  {note 0}
  743.  *
  744.  * When does J(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) take the value of -1,
  745.  * thus satisfy [condition 1]?  The answer depends on P.  Now P is a prime>2,
  746.  * thus P mod 4 == 1 or -1.
  747.  *
  748.  * Take P mod 4 == 1:
  749.  *
  750.  *    P mod 4 == 1  implies  (-1)^((h*2^n-2)*(P-1)/4) == 1
  751.  *
  752.  * Thus:
  753.  *
  754.  *    L(D, h*2^n-1) == L(h*2^n-1 mod P, P) * (-1)^((h*2^n-2)*(P-1)/4)
  755.  *              == L(h*2^n-1 mod P, P)
  756.  *              == J(h*2^n-1 mod P, P)
  757.  *
  758.  * Take P mod 4 == -1:
  759.  *
  760.  *    P mod 4 == -1  implies  (-1)^((h*2^n-2)*(P-1)/4) == -1
  761.  *
  762.  * Thus:
  763.  *
  764.  *    L(D, h*2^n-1) == L(h*2^n-1 mod P, P) * (-1)^((h*2^n-2)*(P-1)/4)
  765.  *              == L(h*2^n-1 mod P, P) * -1
  766.  *              == -J(h*2^n-1 mod P, P)
  767.  *
  768.  * Therefore [condition 1] is met if, and only if, one of the following
  769.  * to cases are true:
  770.  *
  771.  *    P mod 4 ==  1  and  J(h*2^n-1 mod P, P) == -1
  772.  *    P mod 4 == -1  and  J(h*2^n-1 mod P, P) ==  1
  773.  *
  774.  * Now consider [condition 2]:
  775.  *
  776.  *    L(r, h*2^n-1) * (a^2 - b^2*D)/r == -1    [condition 2]
  777.  *
  778.  * We select only a, b, r and D values where:
  779.  *
  780.  *    (a^2 - b^2*D)/r == -1
  781.  *
  782.  * Therefore in order for [condition 2] to be met, we must show that:
  783.  *
  784.  *    L(r, h*2^n-1) == 1
  785.  *
  786.  * If we select r to be of the form:
  787.  *
  788.  *    r == (2^j)*(3^k)*(z^2)            (j>=0, k>=0, z>0)
  789.  *
  790.  * then by use of {note 7}:
  791.  *
  792.  *    L(r, h*2^n-1) == L((2^j)*(3^k)*(z^2), h*2^n-1)
  793.  *              == 1                           {note 2}
  794.  *
  795.  * and thus, [condition 2] is met.
  796.  *
  797.  * If we select r to be of the form:
  798.  *
  799.  *    r == Q*(2^j)*(3^k)*(z^2)        (Q is prime>2, j>=0, k>=0, z>0)
  800.  *
  801.  * then by use of {note 7}:
  802.  *
  803.  *    L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1)
  804.  *              == L(Q, h*2^n-1) * L((2^j)*(3^k)*(z^2), h*2^n-1) {A3.5}
  805.  *              == L(Q, h*2^n-1) * 1                   {note 2}
  806.  *              == L(h*2^n-1, Q) * (-1)^((h*2^n-2)*(Q-1)/4)      {A3.8}
  807.  *              == L(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4)  {note 1}
  808.  *              == J(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4)  {note 0}
  809.  *
  810.  * When does J(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) take the value of 1,
  811.  * thus satisfy [condition 2]?  The answer depends on Q.  Now Q is a prime>2,
  812.  * thus Q mod 4 == 1 or -1.
  813.  *
  814.  * Take Q mod 4 == 1:
  815.  *
  816.  *    Q mod 4 == 1  implies  (-1)^((h*2^n-2)*(Q-1)/4) == 1
  817.  *
  818.  * Thus:
  819.  *
  820.  *    L(D, h*2^n-1) == L(h*2^n-1 mod Q, Q) * (-1)^((h*2^n-2)*(Q-1)/4)
  821.  *              == L(h*2^n-1 mod Q, Q)
  822.  *              == J(h*2^n-1 mod Q, Q)
  823.  *
  824.  * Take Q mod 4 == -1:
  825.  *
  826.  *    Q mod 4 == -1  implies  (-1)^((h*2^n-2)*(Q-1)/4) == -1
  827.  *
  828.  * Thus:
  829.  *
  830.  *    L(D, h*2^n-1) == L(h*2^n-1 mod Q, Q) * (-1)^((h*2^n-2)*(Q-1)/4)
  831.  *              == L(h*2^n-1 mod Q, Q) * -1
  832.  *              == -J(h*2^n-1 mod Q, Q)
  833.  *
  834.  * Therefore [condition 2] is met by selecting  D = Q*(2^j)*(3^k)*(z^2),
  835.  * where Q is prime>2, j>=0, k>=0, z>0; if and only if one of the following
  836.  * to cases are true:
  837.  *
  838.  *    Q mod 4 ==  1  and  J(h*2^n-1 mod Q, Q) == 1
  839.  *    Q mod 4 == -1  and  J(h*2^n-1 mod Q, Q) == -1
  840.  *
  841.  ***
  842.  *
  843.  * In conclusion, we can compute v(1) by attempting to do the following:
  844.  *
  845.  * h mod 3 != 0
  846.  *
  847.  *     we return:
  848.  *
  849.  *       v(1) == 4
  850.  *
  851.  * h mod 3 == 0
  852.  *
  853.  *     define:
  854.  *
  855.  *       r == abs(a^2 - b^2*D)
  856.  *       alpha == (a + b*sqrt(D))^2/r
  857.  *
  858.  *     we return:
  859.  *
  860.  *       v(1) = alpha^1 + alpha^(-1)
  861.  *
  862.  *     if and only if we can find a given a, b, D that obey all the
  863.  *     following selection rules:
  864.  *
  865.  *       D is square free
  866.  *
  867.  *       D == P*(2^f)*(3^g)        (P is prime>2, f,g == 0 or 1)
  868.  *
  869.  *       alpha = epsilon^s        (s>0, epsilon fundamental unit
  870.  *                          of K(sqrt(D)))
  871.  *
  872.  *       r == Q*(2^j)*(3^k)*(z^2)    (Q==1 or Q is prime>2, j>=0, k>=0, z>0)
  873.  *
  874.  *         one of the following is true:
  875.  *           P mod 4 ==  1  and  J(h*2^n-1 mod P, P) == -1
  876.  *           P mod 4 == -1  and  J(h*2^n-1 mod P, P) ==  1
  877.  *
  878.  *       if Q is prime, then one of the following is true:
  879.  *           Q mod 4 ==  1  and  J(h*2^n-1 mod Q, Q) == 1
  880.  *           Q mod 4 == -1  and  J(h*2^n-1 mod Q, Q) == -1
  881.  *
  882.  *     If we cannot find a v(1) quickly enough, then we will give up
  883.  *     testing h*2^n-1.  This does not happen too often, so this hack
  884.  *     is not too bad.
  885.  *
  886.  ***
  887.  *
  888.  * input:
  889.  *    h    h as in h*2^n-1
  890.  *    n    n as in h*2^n-1
  891.  *
  892.  * output:
  893.  *    returns v(1), or -1 is there is no quick way
  894.  */
  895. define
  896. gen_v1(h, n)
  897. {
  898.     local d;    /* the 'D' value to try */
  899.     local val_mod;    /* h*2^n-1 mod 'D' */
  900.     local i;
  901.  
  902.     /*
  903.      * check for case 1
  904.      */
  905.     if (h % 3 != 0) {
  906.         /* v(1) is easy to compute */
  907.         return 4;
  908.     }
  909.  
  910.     /*
  911.      * We will try all 'D' values until we find a proper v(1)
  912.      * or run out of 'D' values.
  913.      */
  914.     for (i=0; i < quickmax; ++i) {
  915.  
  916.         /* grab our 'D' value */
  917.         d = d_qval[i];
  918.  
  919.         /* compute h*2^n-1 mod 'D' quickly */
  920.         val_mod = (h*pmod(2,n%(d-1),d)-1) % d;
  921.  
  922.         /*
  923.          * if 'D' mod 4 == 1, then
  924.          *    (h*2^n-1) mod 'D' can not be a quadratic residue of 'D'
  925.          * else
  926.          *    (h*2^n-1) mod 'D' must be a quadratic residue of 'D'
  927.          */
  928.         if (d%4 == 1) {
  929.             /* D mod 4 == 1, so check for J(D, h*2^n-1) == -1 */
  930.             if (jacobi(val_mod, d) == -1) {
  931.                 /* it worked, return the related v(1) value */
  932.                 return v1_qval[i];
  933.             }
  934.         } else {
  935.             /* D mod 4 == -1, so check for J(D, h*2^n-1) == 1 */
  936.             if (jacobi(val_mod, d) == 1) {
  937.                 /* it worked, return the related v(1) value */
  938.                 return v1_qval[i];
  939.             }
  940.         }
  941.     }
  942.  
  943.     /*
  944.      * This is an example of a more complex proof construction.
  945.      * The code above will not be able to find the v(1) for:
  946.      *
  947.      *            81*2^81-1
  948.      *
  949.      * We will check with:
  950.      *
  951.      *    v(1)=81      D=6557      a=79      b=1      r=31
  952.      *
  953.      * Now, D==79*83 and r=79*2^2.  If we show that:
  954.      *
  955.      *    J(h*2^n-1 mod 79, 79) == -1
  956.      *    J(h*2^n-1 mod 83, 83) == 1
  957.      *
  958.      * then we will satisfy [condition 1].  Observe:
  959.      *
  960.       *    79 mod 4 == -1  implies  (-1)^((h*2^n-2)*(79-1)/4) == -1
  961.       *    83 mod 4 == -1  implies  (-1)^((h*2^n-2)*(83-1)/4) == -1
  962.      *
  963.      *    J(D, h*2^n-1) == J(83, h*2^n-1) * J(79, h*2^n-1)
  964.      *              == J(h*2^n-1, 83) * (-1)^((h*2^n-2)*(83-1)/4) *
  965.      *                 J(h*2^n-1, 79) * (-1)^((h*2^n-2)*(79-1)/4)
  966.      *              == J(h*2^n-1 mod 83, 83) * -1 *
  967.      *                 J(h*2^n-1 mod 79, 79) * -1
  968.      *              ==  1 * -1 *
  969.      *             -1 * -1
  970.      *              == -1
  971.      *
  972.      * We will also satisfy [condition 2].  Observe:
  973.      *
  974.      *    (a^2 - b^2*D)/r == (79^2 - 1^1*6557)/31
  975.      *            == -1
  976.      *
  977.      *    L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1)
  978.      *              == L(79, h*2^n-1) * L(2^2, h*2^n-1)
  979.      *              == L(79, h*2^n-1) * 1
  980.      *              == L(h*2^n-1, 79) * (-1)^((h*2^n-2)*(79-1)/4)
  981.      *              == L(h*2^n-1, 79) * -1
  982.      *              == L(h*2^n-1 mod 79, 79) * -1
  983.      *              == J(h*2^n-1 mod 79, 79) * -1
  984.      *              == -1 * -1
  985.      *              == 1
  986.      */
  987.     if (jacobi( ((h*pmod(2,n%(79-1),79)-1)%79), 79 ) == -1 &&
  988.         jacobi( ((h*pmod(2,n%(83-1),83)-1)%83), 83 ) == 1) {
  989.         /* return the associated v(1)=81 */
  990.         return 81;
  991.     }
  992.  
  993.     /* no quick and dirty v(1), so return -1 */
  994.     return -1;
  995. }
  996.  
  997. /*
  998.  * ldebug - print a debug statement
  999.  *
  1000.  * input:
  1001.  *    funct    name of calling function
  1002.  *    str    string to print
  1003.  */
  1004. define
  1005. ldebug(funct, str)
  1006. {
  1007.     if (dbg == 1) {
  1008.         print "DEBUG:", funct:":", str;
  1009.     }
  1010.     return;
  1011. }
  1012.  
  1013. global lib_debug;
  1014. if (!isnum(lib_debug) || lib_debug>0) print "lucas(h, n) defined";
  1015.